# Load Data
fig.path <- "results/graph/"
fig.name <- paste0(fig.path, "app_fig_d02.pdf")

RegDatMeets <- readRDS("data/regdat_meetings.RDS")

## Regressions

konts <- c( "log_pop" , "SR",
         "AV_AGE_F",   "SC1" , "TFR" ,
         "F_SMAM" , 
         "HC1" ,"HC2" ,  "HC4", "HC5",
         "SINGLE_PER" ,"F_CEL_4554", "M_CEL_4554", "FMAR_PRATE",
         "ECMR","ACRES", "POP_DENS" )

fkont <- paste0("Ftime2treat*", konts)
vars <- c( "Ftime2treat*factor(march_path)", "Ftime2treat", "factor(march_path)")

m1 <-felm(as.formula(RegFor( y = "tot_meet", x = c(vars, fkont, "roman_road") ,
   FE = "county_name" , IV="0", clust = "county_name" )),
   data= RegDatMeets)

summary(m1)
kall <- str_subset(names(coefficients(m1)), pattern = "^Ftime2treat.+TRUE$")
coefs <- coefficients(m1)[kall]
se <- m1$cse[kall]

df <-  data.frame(coefs, se)
names(df) <- c("coef", "se")
df$time2treat <- c(-29:-22, 23:31)
df <- rbind(df, c(0,0,-21))

df$lo10 <- df$coef - 1.65*df$se
df$lo5 <- df$coef - 1.96*df$se
df$up10 <- df$coef + 1.65*df$se
df$up5 <- df$coef + 1.96*df$se

kols <- ghibli_palette(name = "PonyoMedium")[c(2,4)]

p.t <- ggplot(data = df, aes(x = factor(time2treat))) +
geom_pointrange(aes( y= coef, ymin=lo5, ymax = up5) , 
      size = 2, fatten = 1,  alpha = 0.9, position = position_dodge(width = 0.4)) +
geom_linerange(aes( ymin=lo10, ymax = up10) , 
      size = 3,  alpha = 0.35, position = position_dodge(width = 0.4)) +
geom_vline(xintercept= 9,linetype="dashed", size=0.2) +
geom_hline(yintercept=0,linetype="dotted", size=0.2) +
theme_light(base_size = 16) + labs( x = "Weeks to March" , y = "Effect Week x March") +
theme(legend.position = "bottom",
        legend.direction = "vertical") +
guides(shape = guide_legend(override.aes = list(size = 0.5)))

ggsave(fig.name, plot = p.t, height = 6, width = 8 )
rm(list=
    setdiff(ls()[sapply(ls(), function(x) any(class(get(x)) == 'data.frame'))],
    c("Re","Ind")))